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Abstract : 

Using Nonrelativistic QCD on the lattice we studied the mass spectrum of quarkonium 
systems nonperturbatively for a range of the bar quark mass. We determined two products 
of the matrix elements involved in quarkonium decays and studied the mass dependence of 
the results. We predict from our calculations the leptonic decay width of T, and use the 
mass dependence to predict the leptonic decay width of J/tp. These calculations agree with 
the experimental results. In lattice NRQCD an additional parameter n is introduced, and 
we study the sensitivity of our results to the choice of n. 



1. Introduction. In quarkonium the heavy quark moves with a small velocity v 2 , so 
nonrelativistic quantum chromodynamics (NRQCD) may be used as a good approximation 
to describe quarkonium systems. Recent studies have shown that quarkonium systems can 
be well approximated on lattice by solving NRQCD nonperturbatively[l-4]. This opens an 
opportunity for precise tests of QCD. It is important to distinguish these first principles 
calculations from model calculations, for example potential model calculations, which may 
be inspired by QCD, but which depend on phenomenological parameters. 

We study quarkonium systems using NRQCD on a 16 3 x 48 lattice. We use 20 quenched 
gauge configurations with f3 = 6.0. In previous studies [4], the mass spectrum of bottonium 
was extensively studied and shown to agree precisely with experiment [4]. We emphasize 
in our work decay matrix elements, one of which is related to relativistic corrections. 
However in order to set the physical scale, we need to calculate the mass spectrum. We 
have simulated quarkonium systems for a range of the heavy quark mass to study the mass 
dependence of the matrix elements, in the hope that it may be possible to extrapolate 
our results for the matrix elements from bottonium to charmonium. We must approach 
charmonium by extrapolation because the lattice spacing at f3 = 6.0 is too small for direct 
simulation of charmonium with lattice NRQCD [5]. Another aspect of our work is to study 
the n-dependence of our results, where n is a parameter introduced into lattice NRQCD 
for numerical reasons. 

The action of NRQCD for quarks on a lattice can be defined as [6]: 

s Q = ^(x)^ X ) - ^, t + i)(i - ^r^(x,t)(i - §^)*v(x,*)} . (i) 

In Eq.(l) Ho is the Hamitonian operator on the lattice, U^(x) is the gauge link, and n 
is an integer parameter. We only retain the action up to order v 2 . In our work tadpole 
improvement [7] is implemented, i.e., each gauge link U ^ is replaced by U^/uo, where uo = 
(|Tr?7pi a q) 4 . For the configurations used here uq = 0.8778. With this action the quark 
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propagator G(x, t) satisfies the evolution equation 

G(x,f + 1) = k.o^+i.o + (1 - ^) n U}(x,t)(l - f^) n G(x,f) 

An An 

H - A(2) h h - 3(1 (2) 
-Oo — "o, ii-o — 

2M Q Mq 
with G(x,t)=0 for t < 0. In Eq.(2) A^ 2 -* is the lattice Laplacian, and Mq is the quark 
mass parameter in lattice units. We subtracted a constant in Ho } this has an effect that 
the results from the mean field theory for the mass renormalization constant Z m and the 
zero point energy Eo are independent of n and also of Mq . The integer n is introduced to 
avoid the numerical instability when high momentum modes occur [8]. The introduction of 
n has an effect only at order a 2 , where a is the lattice spacing. Therefore one expects that 
results from simulations should not have a strong dependence on n as a — > 0. An estimate[6] 
suggests that n should be larger than at fi = 6.0 to avoid the instability. With this 
estimate n = 1 should be enough large for Mq = 2.0. However we will see that numerical 
instability may still occur for Mq around 2.0 and this has effects in the determination 
of the matrix elements. The numerical values of the matrix elements determined with 
n = 1 and n = 2 are significantly different for Mq < 2.0. For the mass spectrum this 
n-dependence is not significant. 

With the action of Eq.(l) spin-symmetry is an exact symmetry, hence quarkonia with 
the same orbital angular momentum have the same mass. We will only consider S"-wave 
and P-wave quarkonia without radial excitations. 

2. The Mass Spectrum. Since the action in Eq.(l) is nonrelativistic, the absolute 
energy scale is unknown. However the zero-point energy Eo of a quark can be calculated 
perturbatively. The exact mass of a quarkonium state, for example an S"-wave state, is 
related to Eo by: 

M s = 2(a- 1 Z M MQ - Eo) + E s . (3) 
Here Es is the nonrelativistic energy of the quarkonium and it can be measured in lattice 



simulations. We will use Eq.(3) to determine the mass of a S"-wave quarkonium. For Zm 
and Eo we use the results from the mean field theory. With Ho given in Eq.(2) they are: 

E = ciE = - In u , Z M = — (4) 

u 

Although the results for one loop corrections to Eo and Zm hi Eq.(2) exist, we consider it 
is not necessary to include them because we employ NRQCD in Eq.(l) only at the order 
of v 2 and the corrections with the tadpole improvement are small. 
To determine Es we measure the correlator on the lattice: 

H S (t) = Y,(o\xH^tm^Ho)x(om (5) 

X 

In Eq. (5), x( x ) 1S the field for the anti-quark. For large t, Hs(t) takes the asymptotic 
form: 

H s (t)n\(S\xH\0)\ 2 e- tE '. (6) 

From this we determine Es and the matrix element. We measure the correlators with 2 4 
initial points for n = 1 and n = 2. We used the standard fit to extract Es and the matrix 
element. Our results for Ms = aMs are given in table 1. 



Table 1 



Mq 


1.5 


1.7 


2.0 


2.3 


2.6 


M s (n = 1) 


5.2436(3) 


5.5171(5) 


5.9363(5) 


6.5440(3) 


6.9391(4) 


M s (n = 2) 


5.2132(5) 


5.4351(5) 


5.8552(4) 


6.3488(4) 


6.8828(5) 



From Table 1 one can see that the masses determined from n = 1 and n = 2 are 
not significantly different, the differences are only at 1 — 2% in the MQ-range we consider. 
We also measured the correlator for P-wave quarkonium to determine the mass spliting 
between the S"-wave and P-wave quarkonium at Mq = 2.6; the result is 

AMp-s = 0.352(6), for n = 1 

(7) 

AMps = 0.331(8), for n = 2. 
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Given the fact that the mass spliting is not sensitive to the quark mass in the range 
considered here, we use the result with n = 2 in Eq.(7) and use the experimental result for 
bottonium to estimate the lattice spacing. We obtain a -1 = 1.33(3)GeV. In this estimation 
one should keep in mind that the systematic error is larger than the statistical error quoted 
here. One source of the systematic error is the neglect of the terms of order v 4 in the action. 
For bottonium v 2 0.1. This means the systematic error due to the 0(v 4 ) terms is about 
10% for AMp_5, and hence also at least 10% in the estimate of a. With this we conclude 
that the S"-wave quarkonium simulated here at Mq = 2.6 approximately corresponds to 
bottonium with the predicted mass Ms = 9.15(20)GeV. The pole mass of the b-quark can 
also be determined: = 3.9(l)GeV. 

Our results for the spectrum are compatible with the results of previous studies by 
other groups. It is interesting to make a comparson with the results from the precise study 
of [4], where the terms at the order of v 4 are included in the lattice action of NRQCD 
and some improvement reducing the effect of the finite lattice spacing is also made. The 
results from there are that the quarkonium at Mq = 1.71 corresponds to bottonium, the 
pole mass of b quark is 4.7GeV and the inverse of the lattice spacing is 2.4GeV. Comparing 
these with our results above gives a feeling how significant the effect of the accurate action 
of [4] can be in the physical results. The lattice spacing determined here is with a factor of 
2 larger than that of [4]. This is expected since the action employed in [4] is more accurate 
in v than ours and is improved to reduce the effect of the finite lattice spacing. 

3. The Matrix Elements. Recently Bodwin, Braaten and Lepage[9] have treated 
quarkonium systems rigorously within QCD. This contrasts with earlier treatments within 
the potential model. A series of factorized forms for the decay and production rate 
of quarkonia were obtained, where the nonperturbative physics is represented through 
NRQCD matrix elements. With this work one can also systematically account for rela- 
tivistic corrections in decay and production processes. For example, the decay rate for 
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1 So — > 77 and for 3 Si — > £ + £ can be written 



as 



TCSo - 77) = 2 7 rQW{^(0|xWSo)( 1 So|V> t x|0) 

Q 

- ^t^oI^t v) 2 *P\ 1 So)( 1 s \^ x \o)}} 

Q 

rCSo - tn-) = ^^{^(o| x V^| 3 5 1 )( 3 5 1 |V'V a |o) 

Q 

4 Re{(0|xW^ D)V| 3 Si)( 3 Si|V>W|0)}} 



(8) 



3MJ"~ IV " I/V ~' v 2 

at the leading order in a and a s . In Eq.(8) terms with Mq 4 are relativistic corrections 
of order v 2 . The matrix elements in Eq.(8) are defined in NRQCD and they can only 
be calculated nonperturbatively. For bottonium the state 1 So has still not been found 
experimentally. Since the spin-symmetry is an exact symmetry in the approximation used 
here, the two products of the matrix elements in r( 1 S'o — > 77) are equal to those in 
r( 3 So — > £ + £~). They have also a simple relation to the matrix elements in hadronic 
decay. We study these products directly, introducing the notation: 

F 1 = (0\ x m i So)( 1 So\^ X \0) 

-' ~ ( 9 ) 
G 1 = Re{(0| X t(_L nf^So)( 1 So\^ X \0)}. 

The quantity F\ has dimension 3 in mass and G\ has dimension 5 in mass. F\ is propor- 
tional to the square of the absolute value of the wave function at the origin. On the lattice 
F\ can be extracted form the correlation function in Eq.(5). To measure G\ we construct 
a suitable correlation function, in which we use the covariant centered difference on the 
lattice for the covariant derivative D. The quantity G\ is also studied in [10], where the 
same lattice action is used as in Eq.(l). However, it is claimed in [10] that the quarkonium 
simulated at Mq = 1.5 corresponds to bottonium, which is in contrary to our result above. 
Our results for F\ = a 3 F\ and G\ = a 5 G\ are given in Table 2. 

Table 2 
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1 ^ 
1 .o 


1 7 


9 n 


9 ^ 


9 

Z.L) 


9 Q 


F l (n = l) 


2.86(4) 


1.92(4) 


1.58(2) 


1.56(2) 


1.68(2) 


1.88(2) 


F 1 (n = 2) 


1.16(2) 


1.20(1) 


1.35(1) 


1.66(2) 


1.82(2) 


2.10(3) 


G^n = 1) 


1.19(1) 


0.64(1) 


0.45(1) 


0.43(1) 


0.48(1) 


0.58(2) 


G^n = 2) 


0.224(8) 


0.248(6) 


0.31(1) 


0.41(1) 


0.53(1) 


0.67(2) 



We have seen that the mass spectrum is not sensitive to the parameter n, however 
the results for F\ and G\ are quite different for different n. The parameter F\ determined 
with n = 1 decreases linearly with Mq from 2.9 to 2.3, then it increases as Mq decreases 
further. From physical arguments one expects that F\ decreases if Mq decreases. The 
possible reason for the deviation from these expectation can be that re = 1 is still not large 
enough to prevent the numerical instability caused by high momentum modes arround 
Mq = 2.0 in the determination of F\. We conclude that the calculation with n = 1 can 
not give correct results for Mq < 2.0, where local sources are used to calculate quark 
propagators in Eq.(2). In the following we will only take the results with n = 2 for 
discussions. 

From our data F\ is approximately proportional to Mq in the Mq range from 1.7 to 
2.9: F\ ~ 0.7 Mq. This behaviour is expected, since F\ is expected to be proportinal to 
Mq. In our approximation Mq = a -1 Z mMq is proportional to Mq. But at Mq = 1.5, F\ 
deviates from this proportionality relation, as can be seen in Fig.l where the straight line 
represents the relation F\ = 0.7 Mq. Similarly, we find that the values of G\ for Mq from 
2.0 to 2.9 satify the relation G\ = 0.078Mq very well, but the values at Mq = 1.5 and 1.7 
deviate from the relation. We suggest that the reason for the deviations at lower Mq is the 
same as discussed above. In Fig. 2 we plot the results for G\, the curve is G\ = 0.078Mq. 
From Fig.l and Fig. 2 one observes that G\ for Mq from 2.0 to 2.9 fits a parabola much 
better than F\ fits a straight line. This is confirmed by the % 2 of the fits: for G\ x 2 1S 
0.61, while x 2 of the fit for F 1 is 20 for Mq from 1.7, 2.0 and 2.3 to 2.9. If we include the 
data point at Mq = 1.5 into the fit, \ 2 becomes 45! 

6 



With these results one can predict the leptonic decay width of T. We use the following 
formula for doing this: 

16 1 -41 

r( 3 ^ - £+£-) = 2vrQ 2 a 2 {(l - -^(Mg^la" 1 ) 3 ^ - - j^io^ fG, }. (10) 

Here the one loop correction from QCD in the coffecient of F\ term is also included. For T 
we take Mq = = 4.7GeV, a = 1/128 and a s {Mf,) = 0.20, where we used the one-loop 
/3-function and the experimental value for a s (Mz) = 0.115 to evolve a s to /i = We 
obtain T(T — > e + e~) 1.7KeV. Our result is in agreement with the experimental result, 
r(T — > e + e~) 1.3KeV. However, as discussed for the a -1 determination before, a large 
systematic error can occur in our result and it can not be taken as a precise prediction. That 
is why we have not quoted any error in our estimate above. The G\ term gives in Eq.(10) 
a small negative contribution which is only 5% of the decay width. Another way to use 
our lattice results to predict the decay width is to extract directly from lattice simulations 
the dimensionless quantities Fi/Mj = /•', \Z\i\l (} Y' and Gi/Mj = Gi/(Z m Mq) 5 for the 
estimate of the decay width. The advantage is that a direct use of a -1 is avoided. But 
since these quantities depends on a power of the quark mass, a precise location of T in the 
parameter space of Mq is needed for a precise prediction. The decay width determined in 
that way with Mb = 2.6 is 2.8KeV, which is not too far away from either the experimental 
value or our other estimate. 

With the mass dependence of F\ and G\ we can also estimate the leptonic decay width 
of J/V>. Assuming M h : M c = M h : M c and taking M c = 1.3GeV we obtain M c = 0.719. 
For a s at /i = M c we use the experimental value a s (M T ) = 0.355 and evolve this value 
to fi = M c . Using Gi = 0.078MJ and F 1 = 0.7M Q we obtain with Eq.(10): T(J/ip -> 
e + e~) 7.1KeV, which is not too far from the experimental result: T(J/ip — > e + e~) 
5.4KeV. In our estimates the QCD corrections are very important, in particular the factor 
(1 — j^Oi s ) is 0.29 for J/tp, and 0.66 for T, but without the QCD correction the factor is 1. 
We also see that the relativistic correction is large for G\ because of the relatively small c 
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quark mass. This correction is at the level of 20%. Our final result is the determination 
of V(rj c — > 77). With the extrapolated values of F\ and G\ we obtain: T(r] c — > 77) 
23KeV, where the one loop QCD correction is included. Comparing the experimental result 
V(t] c — > 77) 6.6KeV our value is three times too large. The fact that our prediction for 
T(J/ip — > e + e~) is close to the experimental result gives some surport to the validity of 
the mass dependence for F\ and G\ . However, the discrepancy in V(t] c — > 77) indicates 
that the mass dependence may need to be modified in order to get both V(t] c — > 77) 
and T(J/ip — > e + e~) in agreement with experiment. In our work we neglected the 0(v 4 ) 
terms in the action and have taken only the tree-level results for renormalization constants 
to convert lattice results into those of the continuum. One can hope that including the 
0(v 4 ) terms, and by taking higher order effects in renormalization constants into account, 
the results will be improved. The effect of renormalization constants may be significant 
in converting lattice results for F\ to continuum results, since at the one loop level, the 
operator in F\ is mixed with that in G\. All of the relevant renormalization constants are 
not yet available at the one-loop level. 

4. Summary. In this work we have studied the properties of quarkonium systems 
by solving NRQCD on lattice nonperturbatively. The results for the mass spectrum are 
compatible with the results from other groups and with experiment. We studied two 
products of matrix elements involved in quarkonium decays, one of which is related to 
the relativistic corrections. Our prediction of T(T — > e + e~) 1.7KeV is in agreement 
with experiment. It is found that the relativistic correction is small (5%) as expected. It 
should be stressed that this result is a QCD prediction and is not a result from model 
calculations. We simulated quarkniom systems in a range of the quark mass and found 
the mass dependence of F\ and G\. F\ is proportional to Mq as expected and G\ is 
proportional to Mq- With these relations we extrapolate our prediction for T(T — > e + e~) 
to T(J/ip — > e + e - ). The prediction for T(J/ip — > e + e~) is in reasonable agreement with 

8 



experiment. A large relativistic correction (ss 22%) in charmonium is found. However our 
prediction for V(t] c — > 77) is 3 times larger than the experimental value. 

In NRQCD on lattice a extra parameter n is introduced for preventing numerical 
instabilities from high momentum modes. We studied the n-dependence of our results by 
taking n = 1 and n = 2. The mass spectrum does not depend on n significantly. However, 
the difference between the results for F\ and G\ with defferent n is large at Mq < 2.0. In 
addition, the value of F\ determined with re = 1 at Mq < 2.0 increases with decreasing 
quark mass, this is against physical expectations. We think that the results with n = 2 for 
Mq from 2.0 to 2.9 are reasonable. Our prediction for the decay widths is based on these 
results. 
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Figure Caption 

Fig.l The values of F\ vs Mq. The x-axis is for Mq, the y-axis is for F\. The points 
are the data points. 

Fig. 2 The values of G\ vs Mq. The x-axis is for Mq, the y-axis is for G\. The points 
are the data points. 
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